Climbing the Density Functional Ladder: Non-Empirical Meta-Generalized Gradient 
Approximation Designed for Molecules and Solids 
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The electron density, its gradient, and the Kohn-Sham orbital kinetic energy density are the 
local ingredients of a meta-generalized gradient approximation (meta-GGA). We construct a meta- 
GGA density functional for the exchange-correlation energy that satisfies exact constraints without 
empirical parameters. The exchange and correlation terms respect two paradigms: one- or two- 
electron densities and slowly-varying densities, and so describe both molecules and solids with high 
accuracy, as shown by extensive numerical tests. This functional completes the third rung of "Jacob's 
ladder" of approximations, above the local spin density and GGA rungs. 

PACS numbers: 71.15.Mb,31.15.Ew,71.45.Gm 
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Kohn-Sham spin-density functional theory reduces 
the many-electron ground state problem to a self- 
consistent noninteracting-electron form that is exact in 
principle for the density and energy, requiring in prac- 
tice an approximation for the exchange-correlation (xc) 
energy functional E xc [rir, nA. No other method achieves 
comparable accuracy at the same cost. The exact func- 
tional is universal; its approximations should be usefully 
accurate for both molecules and solids, and thus for inter- 
mediate cases (clusters, biological molecules) or combi- 
nations (chemisorption and catalysis on a solid surface). 
The paradigm densities of quantum chemistry (hydro- 
gen atom and electron pair bond) and condensed matter 
physics (uniform electron gas) thus deserve special re- 
spect. Semi-empirical functionals can fail outside their 
fitting sets 0, Q ; those fitted only to molecules can be 
unsuitable for solids. Alternatively, functionals can be 
constructed to satisfy exact constraints on E xc [n-\,n\]. 
Non-empirical functionals are not fitted to actual or com- 
puter experiments for real systems, but are validated by 
such data. 

The first three rungs of "Jacob's ladder" p| of approx- 
imations can be summarized by the formula 



#xcK> n l] = / d r ne xc(n-T'"4' Vn T> v ™b T T> T l)> (!) 



where n(r) = n-|-(r) + nj.(r) is the total density, and 



occup 



^ 2 



(2) 



is the kinetic energy density for the occupied Kohn-Sham 
orbitals ^^(r), which are nonlocal functionals of the den- 
sity n a {v). (We use atomic units where h — m — e 2 = 1.) 
The first rung, found by dropping the Vn CT and r a depen- 
dences in Eq. JTJ, is the local spin density (LSD) approxi- 
mation exact for the uniform electron gas and often 



usefully accurate for solids. The second rung, found by 
dropping only the T a dependence in Eq. (Q, is the gener- 
alized gradient approximation (GGA) |(| 0, (useful for 
molecules as well) . The Perdew-Burke-Ernzerhof (PBE) 
GGA has two non-empirical derivations based on exact 
properties of the xc hole Q and energy 

The third rung of the ladder, the meta-GGA, makes 
use of all the ingredients shown in Eq. (JTJ. The kinetic 
energy density T a (r) recognizes when n a [r) has one- 
electron character by the condition t ct (r) = (r) , where 
r^(r) = \Vn a \ 2 /8n a is the von Weizsacker kinetic en- 
ergy density for real orbitals. Moreover, Eq. can be 
correct to fourth-order in V for a slowly- var ying den- 
sity [ll|. However, prior meta-GGAs 0, El 111! have 
been at least partly empirical, and have not taken full 
advantage of T a (r). The aim of this work is to construct 
a reliable non-empirical meta-GGA to complete the third 
rung of Jacob's ladder. 

Our starting point is the Perdew-Kurth-Zupan-Blaha 
(PKZB) meta-GGA which by design yields the cor- 
rect exchange and correlation energies through second- 
order in V for a slowly-varying density, and the correct 
correlation energy for any one-electron density. PKZB 
has one empirical parameter (fitted to atomization en- 
ergies) in its exchange part, and has demonstrated un- 
expected successes and failures including: (i) an accu- 
rate strong-interaction limit for the correlation energy of 
a spin-unpolarized non-uniform densit y 1121 . (ii) accu- 
rate surface and atomization pi Jill energies, 
(hi) poor equilibrium bond lengths |lj|, andfiv) poor 
description of hydrogen-bonded complexes While 
PKZB correlation requires only minor technical improve- 
ments, PKZB exchange (the root of the bond-length er- 
rors) needs to take into account both paradigm densi- 
ties. Since it is impossible within the meta-GGA form 
to achieve the exact exchange energy for an arbitrary 
one- or two-electron density (i.e., the fully nonlocal self- 
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interaction correction to the Hartree energy), we shall 
instead require that the meta-GGA exchange potential 
be finite at the nucleus for ground-state one- and two- 
electron densities, an exact constraint satisfied by LSD 
but lost in GGA |(|. This is the key idea of our con- 
struction. Because the PKZB meta-GGA is explained in 
detail in Ref. , we focus here on the added features of 
our proposed "TPSS" meta-GGA functional. 

We begin with the exch ang e energy. Because of the 
exact spin-scaling relation [llj i? x [n|,rtj = E x [2n^]/2 + 
E x [2n±]/2, where E x [n) = E x [n/2,n/2], we need con- 
sider only the spin-unpolarized case. Using the uniform 
density-scaling constraint (Eq. (6) of Ref. [llj), the meta- 
GGA (MGGA) can be written as 



E™ GGA [n] 



d 3 r ne™ il (n)F x (p,z), 



(3) 



where e x nlf (n) — — ^(S^n) 1 / 3 is the exchange energy 
per particle of a uniform electron gas. The enhancement 
factor F x (which equals 1 in LSD) is a function of two 
dimensionless inhomogeneity parameters 



p = |Vn| 2 /[4(3^ 2 ) 2 / 3 n 8 / 3 ] 



zero, qb reduces to the dimensionless variable q defined 
in Eq. (12) of Ref. [HI; b will be defined later. 

For a two-electron ground-state density, z = 1 (a = 0) 
and the meta-GGA reduces to GGA form with an en- 
hancement factor F x (s) = F x (p = s 2 ,z = 1). The cor- 
responding exchange potential v x (r) = SE X /Sn(r) has 
a V 2 n term which diverges at a nucleus unless its co- 
efficient, proportional to dF x /ds (Eq. (24) of Ref. |||), 
vanishes there. Since, for a two-electron exponential den- 
sity, s at a nucleus is 0.376, we eliminate this spurious 
divergence (see Fig. 1 of Ref. 0) by requiring that 



dF x /ds 



s=0.376 



= 0. 



(9) 
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All of the above exact constraints imposed on the 
TPSS meta-GGA, beyond those imposed on the PBE 
GGA 0, Ht , concern small dimensionless density deriva- 
tives. We have no reason to modify the large-p behavior 
of the PBE enhancement factor F x ~ 1 + n — n 2 /(fj,p) 
(p — > oo), where fi — 0.21951. Thus we retain this large- 
p limit, which can describe weak binding (lif. 

To satisfy all of the above conditions, we choose x of 
Eq. JHJ to be 



-W _ 



and z = t /t < 1, where r = Y] r a and 
±|Vn| 2 /n. To ensure the Lieb-Oxford bound (Eq. (7) 
of Ref. 0), we choose 



F x 



1 + K — «/(l + x/k), 



(5) 



where k = 0.804 and x(p, z) > will be defined later. (A 
tighter bound 0] would make k = 0.758.) 

For a density that varies slowly over space, F x should 
recover the fourth-order gradient expansion of Svendscn 
and von Barth 16), 



l + H p+ i^- g 2_J3_ 2 6 

8r 2025 y 405^ F v h y ' 



where q = V 2 n/[4(37r 2 ) 2 / 3 n 5 / 3 ] is the reduced Lapla- 
cian. In PKZB, the gradient coefficient D is an empiri- 
cal parameter equal to 0.113, but in TPSS we set D = 
0, the best numerical estimate for this coefficient |l6|. 
The second-order gradient expansion for r (Eq. (4) of 
Ref. 0) makes it possible to recover Eq. (JSJ from a 
Taylor expansion of F x (p, z). Note in particular that 

q b = (9/20)(a - 1)/[1 + ba(a - l)] 1 / 2 + 2p/3, (7) 

an inhomogeneity parameter constructed from p and z, 
tends to the reduced Laplacian q in the slowly-varying 
limit. In Eq. J7J), 



a = (t - 



r-W 



)/ T "»a = (5p/3)(z- 1 - 1) > 0, 



(8) 



where r umf = ^j(37r 2 ) 2 / 3 n 5 / 3 is the uniform-gas kinetic 
energy density. When the parameter b in Eq. Q is set to 
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+2^(|^ +e W 3 |/(l + ^) 2 , (10) 

where the constants c = 1.59096 and e = 1.537 are cho- 
sen to enforce Eq. © and to yield the correct exchange 
energy (—0.3125 hartree) for the exact ground-state den- 
sity of the hydrogen atom. The form of Eq. 1)10(1. which 
is not uniquely constrained, is chosen to make F x smooth 
and monotonic as a function of s for < a < 1. Finally, 
the parameter b = 0.40 in Eq. J7J is chosen to have the 
smallest value that preserves F x as a monotonically in- 
creasing function of s even for fixed a» 1; this is done 
for esthetic reasons, since 6 = produces nearly the same 
results in our molecular tests. 

Figure 1 shows that TPSS has a strong but local ex- 
change enhancement at small s for a = 0. In contrast, 
PKZB has F x w 1 in the small-s region (see the r s = 
curves for a = and 1 in Figs. 1 and 2 of Ref. |2|]), while 
the PBE GGA F x is altogether independent of a. 

We turn now to correlation, making minor refinements 
(independent of TPSS exchange) to PKZB [lj: 
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FIG. 1: TPSS and PBE exchange enhancement factors as 
functions of the reduced gradient s of Eq. The TPSS 

plots (solid lines) are drawn for four different values of a of 
Eq. (JHJ. Ground state two-electron densities have a = 0; 
slowly- varying densities have a ~ 1 and s < 1. 
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The basic ingredient here is e^ BE (nf, nj,, Vn-f , VnjJ |8J 
which already has the correct second-order gradient ex- 
pansion, scales properly under uniform density-scaling to 
the high- and low-density (weak- and strong-interaction) 
limits, is nonpositive, and vanishes as p — > oo. The 
quantity e c could be eJ?' BE (n ., 0, Vn,, 0) as in PKZB, 
but to ensure the exact constraint E c < for all pos- 
sible densities we take e c — max[eJ; BE (ii ff , 0, Vn a , 0), 

BE (ri|, n^, Vri|, Vrij)]. The switch is almost seamless 
in the rare cases where it occurs such as the Li atom 
(see Fig. 2 of Ref. 0), and would probably not oc- 
cur at all for a perfectly-designed GGA. Note also that 
t w = |Vn| 2 /8n was spin- resolved in PKZB. 

For any C*(C,£), £ c MGGA of Eq. (TTJ properly van- 
ishes for a one-electron density (for which r = t w and 
£ = (nj — n\)/n = ±1). Since a self-interaction cor- 
rection should not change the energy of delocalized elec- 
trons, we set C(0,0) = 0.53 (as in PKZB) and d = 2.8 
hartree -1 . This choice leaves the surface correlation en- 
ergy of jellium unchanged from its PBE GGA values over 
the range of valence-electron bulk densities. 

By properly eliminating the self-correlation error for 
spin-unpolarized (£ = 0) densities, the PKZB and 
TPSS approximations yield correct correlation energies 
for atomic densities (with £ = 0) even in the strong- 
interaction limit 01 j where LSD and GGA fail badly. 
In this limit, electrons are kept apart by Coulomb re- 
pulsion, so they "forget" that they are fermions and the 
exchange-correlation energy E xc for a given density n(r) 
becomes independent of the relative spin polarization £ 
(as it almost does in LSD and PBE GGA). To achieve this 



independence over the range < |£| < 0.7 for Gaussian 
one-electron and other densities of uniform £ requires 



C(£, 0) = 0.53 + 0.87C 2 + 0.50C 4 + 2.26C 6 



(13) 



In a monovalent atom like Li, a self-interaction correc- 
tion should zero out the correlation energy density in the 
valence region without creating any additional correla- 
tion energy density in the core-valence overlap region (a 
PKZB error, Fig. 2 of Ref. We achieve this with 



c(C,0 = 



c(C,o) 



{l + £ 2 [(l + C)- 4 / 3 + (W)- 4 / 3 ]/2} 4 



(14) 



where £ = | V£| /2(37r 2 n) 1//3 is large in the overlap region. 
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(£ also appears in Ref. [lij). For a spin-unpolarized den- 
sity, setting d = in Eq. 1)11(1 recovers PKZB. 

Finally, we have made extensive numerical tests of the 
TPSS meta-GGA for atoms, molecules, solids and jellium 
surfaces. Total energies of atoms are exceptionally accu- 
rate, but we focus in Table I on energy differences and 
related properties of greater chemical and physical inter- 
est. In this table, we also report LSD, PBE GGA, PKZB 
meta-GGA, and PBE0 hybrid [2(| results for compari- 
son. Details of our study, and comparisons with other 
functionals, will be presented in later publications 

Except for the jellium surface 0, all our calculations 
are performed self-consistently (although post-PBE re- 
sults would be similar) by the method for r-dependent 
functionals described in Ref. [2lJ using a developmental 
version of the Gaussian program Alternatively, a 

meta-GGA Kohn-Sham potential could be found by the 
optimized effective potential method |2^|. All molecular 
calculations use the large 6-311++G(3df,3pd) Gaussian 
basis set. Performance of the functionals for molecular 
binding properties is assessed by computing atomization 
energies for the 148 molecules (built up from atoms with 
Z < 17) of the G2 test set 24]. For organic molecules 
larger than those in this set, TPSS is much more ac- 
curate [13, El than PKZB. We also studied the ten 
hydrogen-bonded complexes of Ref. ^4Jj reporting the 
errors with respect to M0ller-Plesset (MP2) predictions 
(which agree with experimental values where available), 
and 17 solids (Li, Na, K, Al, C, Si, SiC, BP, NaCl, NaF, 
LiCl, LiF, MgO, Cu, Rh, Pd, Ag in their normal crystal 
structures) calculated with various basis sets. 

Table I shows that the TPSS meta-GGA gives gen- 
erally excellent results for a wide range of systems and 
properties, correcting overestimated PKZB bond lengths 
in molecules, hydrogen-bonded complexes, and ionic 
solids. We attribute this good performance to the sat- 
isfaction of additional exact constraints on i? xc [nf,nj, 
especially the key constraint of Eq. @. Note also that 
PBE0 20] results are comparable to TPSS, but TPSS 
has a practical advantage, since hybrid functionals are 
not easily evaluated for solids. More generally, exact ex- 
change is inaccessible or too slow in many codes. 
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TABLE I: Statistical summary of the errors of four density functionals for various properties of molecules, solids, and surfaces. 
1 kcal/mol = 0.0434 eV = 0.00159 hartree. For jellium, r s = (3/47rn) 1 ' 3 characterizes bulk density. 







Mean value 




Mean absolute 


errors 




Mean error 


Property (units) 


Test set 


of property 


LSD 


PBE 


PBE0 


PKZB 


TPSS 


of TPSS 


Atomization energy EDo (kcal/mol) G2 (148 mols.) 


478 


83.8 


17.1 


5.1 


4.4 


6.2 


5.4 


Ionization potential (eV) 


G2 (86 species) 


10.9 


0.22 


0.22 


0.20 


0.29 


0.23 


-0.11 


Electron affinity (eV) 


G2 (58 species) 


1.4 


0.26 


0.12 


0.17 


0.14 


0.14 


-0.01 


Bond length r e (A) 


96 molecules 


1.56 


0.013 


0.016 


0.010 


0.027 


0.014 


0.014 


Harmonic frequency u e (cm -1 ) 


82 diatomics 


1430 


48.9 


42.0 


43.6 


51.7 


30.4 


-18.7 


H-bond dissoc. energy Do (kcal/ 


mol) 10 complexes 


13.4 


5.8 


1.0 


0.7 


2.9 


0.6 


0.2 


H-bond lengths r e (A) 


11 H-bonds 


2.06 


0.147 


0.043 


0.032 


0.179 


0.021 


0.021 


H-bond angles (deg) 


13 angles 


111 


4.0 


2.6 


1.8 


3.5 


2.0 


2.0 


Lattice constant (A) 


17 solids 


4.34 


0.069 


0.052 




0.078 


0.037 


0.035 


Bulk modulus (GPa) 


17 solids 


124 


14.6 


8.7 




9.2 


9.1 


-1.9 


XC surface energy (erg/cm 2 ) 


r s = 2,4,6 


1245 


22 


55 


39 


5 


13 


-10 



The atomization energy of a molecule can be regarded 
as the extra surface energy of the separated atoms. Meta- 
GGAs are sophisticated enough to reduce the too-large 
LSD atomization energies and increase the too-small LSD 
jellium surface energies, while GGAs reduce both. 

While LSD and PBE GGA are controlled extrapola- 
tions from the slowly-varying limit, TPSS meta-GGA 
is more like a controlled interpolation between slowly- 
varying and one- or two-electron limits. Our functionals 
are nested like Chinese boxes: LSD is inside PBE GGA, 
and PBE GGA is inside TPSS meta-GGA. In a future 
work, we will carry TPSS forward into the fourth rung 
of Jacob's ladder, where an exact-exchange ingredient 
is employed in a global hybrid || like PBE0, a local 
hybrid 26], or a hyper-GGA [4j. A hyper-GGA is exact 
for any one-electron density, and for any density scaled 
uniformly to the high-density limit. 
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